from osgeo import gdal,ogr,osr,gdalconst
from multiprocessing.pool import ThreadPool
from threading import Thread
from tqdm import tqdm
import queue
import copy
import zipfile
import shutil
import glob
import time
import os


# 使用 UnaryUnion 计算并集
# union_geometry = ogr.Geometry(ogr.wkbGeometryCollection)
# for geom in geometries:
#     union_geometry.AddGeometry(geom)
# union_geometry.UnaryUnion()

# layer.CommitTransaction()  
# dataset_mem = driver.CopyDataSource(dataset, "")
# current_union = geometry.Clone()
# layer.ResetReading()
# layer.CreateSpatialIndex()
# index = layer.GetSpatialIndex()
# for i in range(layer_src.GetFeatureCount()):layer_rst.DeleteFeature(i)
# layer_src.Clip(layer_mtd,layer_rst,callback=gdal.TermProgress_nocb)
# layer.SetAttributeFilter(f"class='{value}'")
# attributes = feature.items()
# layer.SyncToDisk()
# wkt_geometry = geometry.ExportToWkt()
# for field in src_layer.schema:dst_layer.CreateField(field)
# geometry = ogr.CreateGeometryFromWkt(wkt)
#  driver = ogr.GetDriverByName('Memory') # 'ESRI Shapefile'
# mem_ds = driver.CreateDataSource('')
# mem_lyr = mem_ds.CreateLayer('layer',target_proj,geom_type=ogr.wkbPolygon)


def repair_geometry(
    vector_source,
    vector_result,
    layer_name,
):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8 ISO-8859-1
    if isinstance(vector_source,str):
        dataset_src = ogr.Open(vector_source,gdal.GA_ReadOnly)
    else:
        dataset_src = vector_source
    layer_src = dataset_src.GetLayer()
    if isinstance(vector_result,str):
        driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
        dataset_rst = driver.CreateDataSource(vector_result)
    else:
        dataset_rst = vector_result
    layer_rst = dataset_rst.CopyLayer(layer_src,layer_name)
    dataset_src.Destroy()
    with tqdm(total=layer_rst.GetFeatureCount()) as bar:
        for i,feature in enumerate(layer_rst):
            # print(feature.GetField('qhmc'))
            while not(geometry.IsSimple() 
                and geometry.IsValid() 
                and not geometry.IsEmpty()
                ):
                geometry = geometry.MakeValid()
                geometry = geometry.Normalize()
                geometry = geometry.Simplify(0.1)
                geometry = geometry.Buffer(0)
                geometry = geometry.Union(geometry)
            feature.SetGeometry(geometry)
            layer_rst.SetFeature(feature)
            bar.update(1)
    # layer_rst.SyncToDisk()
    # dataset_rst.Destroy()         
    return dataset_rst


# Spatial coordinate reference system
def transform_spatialref(
    vector_source,
    vector_result,
    layer_name,
    projection=None,
    epsg_code=None,
):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8 ISO-8859-1
    if isinstance(vector_source,str):
        dataset_src = ogr.Open(vector_source,gdal.GA_ReadOnly)
    else:
        dataset_src = vector_source
    layer_src = dataset_src.GetLayer()
    if isinstance(vector_result,str):
        driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
        dataset_rst = driver.CreateDataSource(vector_result)
    else:
        dataset_rst = vector_result
    layer_src = dataset_src.GetLayer()
    proj = layer_src.GetSpatialRef()
    if epsg_code is not None:
        projection = osr.SpatialReference()
        projection.ImportFromEPSG(int(epsg_code))  # 指定目标EPSG代码
    transform = osr.CoordinateTransformation(
        layer_src.GetSpatialRef(), projection)
    layer_rst = dataset_rst.CreateLayer(
        layer_name,projection,geom_type=ogr.wkbPolygon)
    for field in layer_src.schema:
        layer_rst.CreateField(field)
    layerdefn = layer_rst.GetLayerDefn()
    for feature in tqdm(layer_src,desc="坐标系转换"):
        geometry = feature.GetGeometryRef()
        geometry.Transform(transform)  # 执行坐标转换
        new_feature = ogr.Feature(layerdefn)
        for i in range(feature.GetFieldCount()):
            new_feature.SetField(i, feature.GetField(i))
        new_feature.SetGeometry(geometry)
        layer_rst.CreateFeature(new_feature)
    # layer_rst.SyncToDisk()
    dataset_src.Destroy()    
    # dataset_rst.Destroy()   
    return dataset_rst


def single_thread(function,static_parameters,task_queue):
    data = []
    while task_queue.empty():
        element = task_q.get()
        result = funtion(*static_parameters,*element)
        data.append(result)
    return data


def is_arithmetic_sequence(lst):
    if len(lst) <= 2:
        return True
    differences = [lst[i + 1] - lst[i] for i in range(len(lst) - 1)]
    return all(d == differences[0] for d in differences)


def check_geometry(shp_path):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    dataset = ogr.Open(shp_path,gdal.GA_ReadOnly)
    layer = dataset.GetLayer()
    layer.ResetReading()
    for feature in tqdm(layer,desc="修复Geometry",disable=False):
        if feature.GetFID() == 12059:
            geometry = feature.GetGeometryRef()
            # if not(geometry.IsSimple() 
            #     and geometry.IsValid() 
            #     and not geometry.IsEmpty()):
            #     print(feature.GetFID())
            print(geometry.GetGeometryName())
            if geometry.GetGeometryName() == "POLYGON":
                points = [point for ring in geometry for point in ring]
            elif geometry.GetGeometryName() == "MULTIPOLYGON":
                points = [point for polygon in geometry for ring in polygon for point in ring]
            print(points)
            # if len(set(tuple(p) for p in points)) < len(points):
            #     print(feature.GetFID())
    return


class Statistics_Zzjg:
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8 ISO-8859-1
    driver = ogr.GetDriverByName('Memory') # 'ESRI Shapefile'
    task_queue = queue.Queue()
    pool_size = 15
    min_area = 5
    field_dict = {
        "nf": [ogr.OFTString,4],
        "qhdm": [ogr.OFTString,20],
        "qhmc":[ ogr.OFTString,20],
        "zwdm": [ogr.OFTString,2],
        "zwmc": [ogr.OFTString,10],
        "mj": [ogr.OFTReal, 24],
        "jcsj": [ogr.OFTDate,None],
        "rspdid": [ogr.OFTString,50]}
    data = []
    
    def __init__(self,jy_shp_path,qh_shp_path,save_dir,layer_name,kwargs):
        dataset = ogr.Open(jy_shp_path,gdal.GA_ReadOnly)
        self.dataset_jy = self.driver.CopyDataSource(dataset, "")
        dataset.Destroy()
        dataset = ogr.Open(qh_shp_path,gdal.GA_ReadOnly)
        self.dataset_qh = self.driver.CopyDataSource(dataset, "")
        dataset.Destroy()
        self.save_dir = save_dir
        self.layer_name = layer_name
        self.kwargs = kwargs

    @staticmethod
    def split_layer(dataset,layercount):
        layer = dataset.GetLayer()
        featurecount = layer.GetFeatureCount()
        step = int((featurecount%layercount+featurecount)/layercount)
        args = [[i,i+step] for i in range(0,featurecount,step)]
        for index,(start_fid,end_fid) in enumerate(args):
            layer.ResetReading()
            layer.SetAttributeFilter(f"{start_fid}<=FID and FID<{end_fid}")
            dataset.CopyLayer(layer,f"{index+1}")
        del layer
        return
    
    @staticmethod
    def repair_geometry(dataset,layer_name):
        layer = dataset.GetLayerByName(layer_name)
        for feature in tqdm(layer,desc="修复Geometry"):
            geometry = feature.GetGeometryRef()
            while not(geometry.IsSimple() 
                and geometry.IsValid() 
                and not geometry.IsEmpty()):
                geometry = geometry.MakeValid()
                geometry = geometry.Normalize()
                geometry = geometry.Union(geometry)
                # geometry = geometry.Simplify(0.1)
                # geometry = geometry.Buffer(0)
            feature.SetGeometry(geometry)
            layer.SetFeature(feature)
        layer.CommitTransaction()
        return 
    
    @classmethod
    def multi_thread_repair_geometry(cls,dataset,layercount):
        threads = []
        for index in range(1,layercount+1): 
            threads.append(Thread(
                name=f"修复geometry:{index}",
                target=cls.repair_geometry,
                args=(dataset,f"{index}")))
        for thread in threads:thread.start()
        for thread in threads:thread.join()
        return
    
    @staticmethod
    def replace_feature(dataset,layercount):
        layer = dataset.GetLayer(0)
        layer.ResetReading()
        layer.SetAttributeFilter(f"0<=FID") 
        for index in range(1,layercount+1):
            layer_temp = dataset.GetLayerByName(f"{index}")
            for feature in layer_temp:
                layer.SetFeature(feature)
            dataset.DeleteLayer(f"{index}")
        return
    
    @classmethod
    def repair_layer(cls,dataset,layercount):
        layer = dataset.GetLayer()
        featurecount = layer.GetFeatureCount()
        if featurecount <= 100000: layercount = 1
        cls.split_layer(dataset,layercount)
        cls.multi_thread_repair_geometry(dataset,layercount)
        cls.replace_feature(dataset,layercount)
        return 

    @staticmethod
    def srstransformer(dataset_src,dataset_mtd):
        layer_src = dataset_src.GetLayer()
        layer_mtd = dataset_mtd.GetLayer()
        transform = osr.CoordinateTransformation(
            layer_src.GetSpatialRef(), 
            layer_mtd.GetSpatialRef())
        for feature in tqdm(layer_src,desc="转换SpatialRef"):
            geometry = feature.GetGeometryRef()
            geometry.Transform(transform)  # 执行坐标转换
            feature.SetGeometry(geometry)
            layer_src.SetFeature(feature)
        layer_src.CommitTransaction()
        return
    
    def progress_bar(self):
        histroy = 0
        total = self.task_queue.qsize()
        with tqdm(total=total,desc="任务进度") as bar:
            while not self.task_queue.empty():
                temp = total- self.task_queue.qsize()
                refresh = temp - histroy
                bar.update(refresh)
                histroy = temp
        return
    
    def statistical_area(self,args):
        driver = ogr.GetDriverByName('Memory') # 'ESRI Shapefile'
        dataset = driver.CopyDataSource(self.dataset_jy,"")
        layer = dataset.GetLayer()
        data = []
        while not self.task_queue.empty():
            qhmc, qhdm, geometry_mtd = self.task_queue.get()
            layer.ResetReading()
            layer.SetSpatialFilter(geometry_mtd)
            area =  0
            for feature in layer:
                geometry = feature.GetGeometryRef()
                intersection = geometry.Intersection(geometry_mtd)
                area = area + intersection.GetArea() * 0.0015
            # if area >= self.min_area :
            data.append({"qhmc":qhmc, "qhdm":qhdm, "mj":area})
        dataset.Destroy()
        return data
    
    def save_data(self):
        gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
        gdal.SetConfigOption('SHAPE_ENCODING', 'GBK') # 解决 SHAPE 文件的属性值 GBK UTF-8
        driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
        dataset = driver.CreateDataSource(self.save_dir)
        layer = dataset.CreateLayer(self.layer_name)
        for field_name,values in self.field_dict.items():
            value_type, value_limit = values
            field_definition = ogr.FieldDefn(field_name, value_type)
            if value_type == ogr.OFTString:
                field_definition.SetWidth(value_limit)
            elif value_type == ogr.OFTReal:
                field_definition.SetPrecision(value_limit)
            # field_definition.SetDefault(nf)
            layer.CreateField(field_definition)
        layerdefn = layer.GetLayerDefn()
        for item in tqdm(self.data,desc="导出数据到dbf文件"):
            feature = ogr.Feature(layerdefn)
            for field_name,field_value in self.kwargs.items():
                feature.SetField(field_name,field_value)
            for field_name,field_value in item.items():
                feature.SetField(field_name,field_value)
            layer.CreateFeature(feature)
        layer.SyncToDisk()
        dataset.Destroy()  
        shp_path = os.path.join(self.save_dir,f"{self.layer_name}.shp")
        shx_path = os.path.join(self.save_dir,f"{self.layer_name}.shx")
        prj_path = os.path.join(self.save_dir,f"{self.layer_name}.prj")
        dbf_path = os.path.join(self.save_dir,f"{self.layer_name}.dbf")
        if os.path.exists(shp_path):os.remove(shp_path)
        if os.path.exists(shx_path):os.remove(shx_path)
        if os.path.exists(prj_path):os.remove(prj_path)
        return dbf_path
    
    def running(self,layercount=8):
        thread_jy = Thread(
            name="解译geometry修复",
            target=self.repair_layer,
            args=(self.dataset_jy,layercount))
        thread_jy.start()
        thread_jy.join()
        
        # thread_qh = Thread(
        #     name="区划geometry修复",
        #     target=self.repair_layer,
        #     args=(self.dataset_qh,layercount))
        # thread_qh.start()
        # thread_qh.join()
        
        # thread_qh = Thread(
        #     name="区划spatialref转换",
        #     target=self.srstransformer,
        #     args=(self.dataset_qh,self.dataset_jy))
        # thread_qh.start()
        # thread_qh.join()
        
        driver = ogr.GetDriverByName('ESRI Shapefile') 
        driver.CopyDataSource(self.dataset_jy, r"/data/fengyy/dateset/temp/0429/2024江西早稻")
        # driver.CopyDataSource(self.dataset_qh, r"/data/fengyy/dateset/temp/0429/2024江西早稻")
        
        # layer_qh = self.dataset_qh.GetLayer()
        # layer_qh.ResetReading()
        # for feature in tqdm(layer_qh,desc="载入数据"):
        #     qhmc = feature.GetField('qhmc')
        #     qhdm = feature.GetField('qhdm')
        #     fid = feature.GetFID()
        #     geometry = feature.GetGeometryRef()
        #     geometry = copy.deepcopy(geometry)
        #     self.task_queue.put([qhmc,qhdm,geometry])
        # Thread(name="进度条线程",target=self.progress_bar,daemon=True).start()
        # iterable = [_ for _ in range(self.pool_size)]
        # pool = ThreadPool(processes=self.pool_size)    
        # data = pool.map(self.statistical_area,iterable=iterable,chunksize=1)
        # pool.close()
        # pool.join()
        # self.dataset_qh.Destroy()
        # self.dataset_jy.Destroy()
        # self.data = [j for i in data for j in i if j['mj'] > self.min_area ]
        # dbf_path = self.save_data()
        # return dbf_path


def shp_clip_tif(dsttif,srctif,shppath):
    layername = os.path.splitext(os.path.basename(shppath))[0]
    # 创建裁剪后的栅格文件
    gdal.Warp(
        dsttif, 
        srctif, 
        format="GTiff", # JPEG
        cutlineDSName=shppath,
        cutlineLayer=layername,
        cutlineWhere=f"countycode='{411521}'",  # 根据要素的属性进行裁剪，这里假设要裁剪的要素的FID为1
        # cutlineBlend=True,
        # cropToCutline=True,
        resampleAlg=gdal.GRIORA_NearestNeighbour, # GRIORA_NearestNeighbour,GRIORA_Bilinear
        # dstNodata=0,
        warpMemoryLimit=1024*1024*1024*50,
        multithread=10,
        callback=gdal.TermProgress_nocb)
    return


def transform_srs_zipfile(src_shp_path,epsg_code,save_dir,layer_name):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'GBK') # 解决 SHAPE 文件的属性值
    dst_shp_path = os.path.join(save_dir,layer_name)
    src_dataset = ogr.Open(src_shp_path)
    src_layer = src_dataset.GetLayer()
    # layer_name = src_layer.GetName()
    source_proj = src_layer.GetSpatialRef()
    target_proj = osr.SpatialReference()
    target_proj.ImportFromEPSG(int(epsg_code))  # 指定目标EPSG代码
    transform = osr.CoordinateTransformation(source_proj, target_proj)
    driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
    dst_dataset = driver.CreateDataSource(dst_shp_path)
    dst_layer = dst_dataset.CreateLayer(
        layer_name,target_proj,geom_type=ogr.wkbPolygon)
    for field in src_layer.schema:dst_layer.CreateField(field)
    dst_layer_defn = dst_layer.GetLayerDefn()
    for feature in tqdm(src_layer,desc="坐标系转换"):
        geometry = feature.GetGeometryRef()
        geometry.Transform(transform)  # 执行坐标转换
        new_feature = ogr.Feature(dst_layer_defn)
        for i in range(feature.GetFieldCount()):
            new_feature.SetField(i, feature.GetField(i))
        new_feature.SetGeometry(geometry)
        dst_layer.CreateFeature(new_feature)
    src_dataset.Destroy() 
    dst_dataset.Destroy()
    paths = glob.glob(f"{dst_shp_path}/*")
    zip_path = dst_shp_path+".zip"
    with zipfile.ZipFile(zip_path, 'w' ,zipfile.ZIP_DEFLATED) as zipf:
        for file_path in paths:
            arcname = os.path.relpath(file_path, save_dir)
            zipf.write(file_path,arcname=arcname)
    shutil.rmtree(dst_shp_path)
    return zip_path


def create_shp(vector_dir,layer_name,data,nf,zwdm,zwmc,jcsj,rspdid):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataset = driver.CreateDataSource(vector_dir)
    layer = dataset.CreateLayer(layer_name)
    field_definition = ogr.FieldDefn("nf", ogr.OFTString)
    field_definition.SetWidth(4)
    field_definition.SetDefault(nf)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("qhdm", ogr.OFTString)
    field_definition.SetWidth(20)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("qhmc", ogr.OFTString)
    field_definition.SetWidth(100)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("zwdm", ogr.OFTString)
    field_definition.SetWidth(2)
    field_definition.SetDefault(zwdm)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("zwmc", ogr.OFTString)
    field_definition.SetWidth(100)
    field_definition.SetDefault(zwmc)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("mj", ogr.OFTReal)
    field_definition.SetPrecision(24)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("jcsj", ogr.OFTDate)
    field_definition.SetDefault(rspdid)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("rspdid", ogr.OFTString)
    field_definition.SetWidth(100)
    field_definition.SetDefault(rspdid)
    layer.CreateField(field_definition)
    layerdefn = layer.GetLayerDefn()
    for item in tqdm(data,desc="导出数据到dbf文件"):
        feature = ogr.Feature(layerdefn)
        feature.SetField("qhdm",item['qhdm'])
        feature.SetField("qhmc",item['qhmc'])
        feature.SetField("mj",item['mj'])
        feature.SetField("nf",nf)
        feature.SetField("zwmc",zwmc)
        feature.SetField("zwdm",zwdm)
        feature.SetField("jcsj",jcsj)
        feature.SetField("rspdid",rspdid)
        layer.CreateFeature(feature)
    dataset.Destroy()   
    return 


def geom_intersection_layer_2_area(args):
    wktjson,geometry = args
    # geometry = ogr.CreateGeometryFromWkt(wkt)
    # dataset = ogr.Open(shp_path)
    layer = dataset.GetLayer()
    
    driver = ogr.GetDriverByName('Memory')
    dataset = driver.CreateDataSource('')
    
    # dataset.CopyLayer(layer,"temp")
    # src_dataset.Destroy()
    # layer = dataset.GetLayer(0)
    layer.SetSpatialFilter(geometry)
    area = 0
    if layer.GetFeatureCount() >= 0:
        for feature in layer:
            geometry = feature.geometry()
            # print(geometry)
            # if not(geometry.IsSimple() 
            #     and geometry.IsValid() 
            #     and not geometry.IsEmpty()
            # ):
                # print(geometry)
            #     geometry = geometry.MakeValid()
            #     geometry = geometry.Normalize()
            # intersection = geometry.Intersection(geometry)
            # area = area + intersection.GetArea() * 0.0015
    layer.ResetReading()
    dataset.Destroy()
    # del layer
    # return area


def loading_vector_dataset(vector_path):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    dataset = ogr.Open(vector_path,gdal.GA_ReadOnly)
    layer = dataset.GetLayer()
    driver = ogr.GetDriverByName('Memory') # 'ESRI Shapefile'
    dataset_mem = driver.CreateDataSource('')
    layer_mem = dataset.CopyLayer(layer,"mem")
    dataset.Destroy()
    # layer_mem.ResetReading()
    # layer_mem.SetSpatialFilter(geometry_mtd)
    return dataset_mem


def stat_cyqh_zzjg(dataset_mem,qhmc,qhdm,geometry_mtd):
    layer_mem.ResetReading()
    layer_mem.SetSpatialFilter(geometry_mtd)
    area = 0
    for feature in layer_mem:
        geometry = feature.GetGeometryRef()
        intersection = geometry.Intersection(geometry_mtd)
        area = area + intersection.GetArea() * 0.0015
    data = {"qhmc":qhmc,"qhdm":qhdm,"mj":area}
    dataset_mem.Destroy()
    return data


def progress_bar(task_queue,desc="任务进度"):
    histroy = 0
    total = task_queue.qsize()
    with tqdm(total=total,desc=desc) as bar:
        while not task_queue.empty():
            temp = total- task_queue.qsize()
            refresh = temp - histroy
            bar.update(refresh)
            histroy = temp
    return


def worker(function,static_parameters,task_queue):
    data = []
    
    
    while task_queue.empty():
        element = task_q.get()
        result = funtion(*static_parameters,*element)
        data.append(result)
    return data


def high_concurrency(
    function,
    static_parameters,
    dynamic_parameters,
    pool_size=15,
    chunksize=1,
    desc="任务进度",
):  
    task_queue = queue.Queue()
    for parameters in dynamic_parameters:
        task_queue.put(parameters)
    thread = threading.Thread(
        target=progress_bar,args=(task_queue,desc),daemon=True)
    thread.start()
    parameters = [function,static_parameters,task_queue] * pool_size
    pool = ThreadPool(pool_size)    
    data = pool.map(worker,parameters,chunksize)
    pool.close()
    pool.join()
    data = [j for i in data for j in i]
    return data


def cy_stat_jy(
    vector_source_path,
    vector_method_path,
    vector_dir,
    layer_name,
    params,
    min_mj=5,
    poolsize=15,
):
    task_q = queue.Queue()
    dataset = ogr.Open(vector_method_path,gdal.GA_ReadOnly)
    layer= dataset.GetLayer()
    count = layer.GetFeatureCount()
    for feature in tqdm(layer):
        qhmc = feature.GetField('qhmc')
        qhdm = feature.GetField('qhdm')
        geometry = feature.GetGeometryRef()
        geometry = copy.deepcopy(geometry)
        task_q.put([qhmc,qhdm,geometry])
    dataset.Destroy()
    
    args = [[vector_source_path,task_q]] * poolsize
    
    thread = threading.Thread(
        target=progress_bar,
        args=(task_q,count),
        daemon=True
    )
    thread.start()
    pool = ThreadPool(poolsize)    
    data = pool.map(cy_stat_jy_sub, args)
    pool.close()
    pool.join()
    
    nf,zwdm,zwmc,jcsj,rspdid = params
    data = [j for i in data for j in i if j["mj"]>min_mj]
    dbf_path = create_shp(
        vector_dir,layer_name,data,nf,zwdm,zwmc,jcsj,rspdid)
    return


def shp_clip_shp(vector_source,vector_method):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8
    if isinstance(vector_source,str):
        dataset_src = ogr.Open(vector_source,gdal.GA_ReadOnly)
    else:
        dataset_src = vector_source
    if isinstance(vector_method,str):
        dataset_mtd = ogr.Open(vector_method,gdal.GA_ReadOnly)
    else:
        dataset_mtd = vector_method
    layer_src = dataset_src.GetLayer()
    layer_mtd = dataset_mtd.GetLayer()
    data = []
    for feature_mtd in tqdm(layer_mtd):
        geometry_mtd = feature_mtd.GetGeometryRef()
        layer_src.ResetReading()
        layer_src.SetSpatialFilter(geometry_mtd)
        qhmc = feature_mtd.GetField('qhmc')
        qhdm = feature_mtd.GetField('qhdm')
        area = 0
        for feature_src in layer_src:
            geometry_src = feature_src.GetGeometryRef()
            intersection = geometry_src.Intersection(geometry_mtd)
            area = area + intersection.GetArea() * 0.0015
        if layer_src.GetFeatureCount()!=0 and area>= 5: 
            data.append([qhmc,qhdm,area])
    dataset_src.Destroy()
    dataset_mtd.Destroy()
    return data


def layer_clip_layer(
    vector_source,
    vector_method,
    vector_result,
    layer_name,
    filter_params=None
):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8
    if isinstance(vector_source,str):
        dataset_src = ogr.Open(vector_source,gdal.GA_ReadOnly)
    else:
        dataset_src = vector_source
    if isinstance(vector_method,str):
        dataset_mtd = ogr.Open(vector_method,gdal.GA_ReadOnly)
    else:
        dataset_mtd = vector_method
    if isinstance(vector_result,str):
        driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
        dataset_rst = driver.CreateDataSource(vector_result)
    else:
        dataset_rst = vector_result
    layer_src = dataset_src.GetLayer()
    layer_mtd = dataset_mtd.GetLayer()
    if filter_params is not None:
        layer_mtd.SetAttributeFilter(filter_params)
    with tqdm(total=layer_mtd.GetFeatureCount()) as bar:
        for i,feature in enumerate(layer_mtd):
            geometry = feature.GetGeometryRef()
            if i == 0:
                current_union = geometry.Clone()
            current_union = current_union.Union(geometry).Clone()
            bar.update(1)
    layer_src.SetSpatialFilter(current_union)
    layer_rst = dataset_rst.CopyLayer(layer_src,layer_name)
    layer_rst.SyncToDisk()
    dataset_src.Destroy()
    dataset_mtd.Destroy()
    dataset_rst.Destroy()
    return 


def stat_field(
    vector_source,
    field_name
):  
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8 ISO-8859-1
    if isinstance(vector_source,str):
        dataset = ogr.Open(vector_source,gdal.GA_ReadOnly)
    else:
        dataset = vector_source
    layer = dataset.GetLayer()
    data = []
    layer.ResetReading()
    for feature in tqdm(layer):
        value = feature.GetField(field_name)
        data.append(value)
    dataset.Destroy()   
    return set(data)


def fusion_geometry(
    vector_source,
    filter_params,
    vector_result,
    layer_name,
):  
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'UTF-8') # 解决 SHAPE 文件的属性值 GBK UTF-8 ISO-8859-1
    if isinstance(vector_source,str):
        dataset_src = ogr.Open(vector_source,gdal.GA_ReadOnly)
    else:
        dataset_src = vector_source
    layer_src = dataset_src.GetLayer()
    spatialRef = layer_src.GetSpatialRef()
    layer_src.SetAttributeFilter(filter_params)
    if isinstance(vector_result,str):
        driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
        dataset_rst = driver.CreateDataSource(vector_result)
    else:
        dataset_rst = vector_result
    layer_rst = dataset_rst.CreateLayer(
        layer_name,spatialRef,geom_type=ogr.wkbPolygon)
    layerdefn = layer_rst.GetLayerDefn()
    with tqdm(total=layer_src.GetFeatureCount()) as bar:
        for i,feature in enumerate(layer_src):
            geometry = feature.GetGeometryRef()
            geometry = geometry.Buffer(0)
            if i == 0:
                union_geometry = geometry.Clone()
            else:
                union_geometry = union_geometry.Buffer(0)
                union_geometry = union_geometry.Union(geometry)
                while not(union_geometry.IsSimple() 
                    and union_geometry.IsValid() 
                    and not union_geometry.IsEmpty()
                ):
                    union_geometry = union_geometry.MakeValid()
                    union_geometry = union_geometry.Normalize()
                    union_geometry = union_geometry.Simplify(0.1)
                    union_geometry = union_geometry.Buffer(0)
                    union_geometry = union_geometry.Union(union_geometry)
            bar.update(1)       
    feature = ogr.Feature(layerdefn)
    feature.SetGeometry(union_geometry)
    layer_rst.CreateFeature(feature)
    layer_rst.SyncToDisk()
    dataset_src.Destroy()  
    dataset_rst.Destroy()  
    return



def data_2_dbf(data,savedir,layername,kwargs):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'GBK') # 解决 SHAPE 文件的属性值 GBK UTF-8
    driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
    dataset = driver.CreateDataSource(savedir)
    layer = dataset.CreateLayer(layername)
    field_definition = ogr.FieldDefn("nf", ogr.OFTString)
    field_definition.SetWidth(4)
    field_definition.SetDefault(nf)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("qhdm", ogr.OFTString)
    field_definition.SetWidth(20)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("qhmc", ogr.OFTString)
    field_definition.SetWidth(100)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("zwdm", ogr.OFTString)
    field_definition.SetWidth(2)
    field_definition.SetDefault(zwdm)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("zwmc", ogr.OFTString)
    field_definition.SetWidth(100)
    field_definition.SetDefault(zwmc)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("mj", ogr.OFTReal)
    field_definition.SetPrecision(24)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("jcsj", ogr.OFTDate)
    field_definition.SetDefault(rspdid)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn("rspdid", ogr.OFTString)
    field_definition.SetWidth(100)
    field_definition.SetDefault(rspdid)
    layer.CreateField(field_definition)
    layerdefn = layer.GetLayerDefn()
    for qhmc,qhdm,area in tqdm(data,desc="导出数据到dbf文件"):
        feature = ogr.Feature(layerdefn)
        feature.SetField("qhdm",qhdm)
        feature.SetField("qhmc",qhmc)
        feature.SetField("mj",area)
        feature.SetField("nf",kwargs['nf'])
        feature.SetField("zwmc",kwargs['zwmc'])
        feature.SetField("zwdm",kwargs['zwdm'])
        feature.SetField("jcsj",kwargs['jcsj'])
        feature.SetField("rspdid",kwargs['rspdid'])
        layer.CreateFeature(feature)
    layer.SyncToDisk()
    dataset.Destroy()  
    shp_path = os.path.join(savedir,f"{layername}.shp")
    shx_path = os.path.join(savedir,f"{layername}.shx")
    prj_path = os.path.join(savedir,f"{layername}.prj")
    dbf_path = os.path.join(savedir,f"{layername}.dbf")
    if os.path.exists(shp_path):os.remove(shp_path)
    if os.path.exists(shx_path):os.remove(shx_path)
    if os.path.exists(prj_path):os.remove(prj_path)
    return dbf_path


def zwcy_2_tb(
    src_shp_path,
    dst_shp_path,
    res_shp_path,
    nf,
    zwdm,
    zwmc,
    jcsj,
    rspdid
):  
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'GBK') # 解决 SHAPE 文件的属性值
    
    dataset_jy = ogr.Open(src_shp_path)
    layer_jy = dataset_jy.GetLayer()
    target_proj = layer_jy.GetSpatialRef()
   
    dataset_jy = repair_geometry(dataset_jy)
    
    dataset_qh = transform_srs(dst_shp_path,target_proj)
    dataset_qh = repair_geometry(dataset_qh)
    
    layer_jy = dataset_jy.GetLayer()
    layer_qh = dataset_qh.GetLayer()
    
    driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
    dataset_jy_re = driver.CreateDataSource("/data/fengyy/dateset/temp")
    dataset_jy_re.CopyLayer(layer_jy,"江西水稻")
    
    dataset_jy_re = driver.CreateDataSource("/data/fengyy/dateset/temp")
    dataset_jy_re.CopyLayer(layer_qh,"江西村域")
    
    # print(layer_jy.GetFeatureCount(),layer_qh.GetFeatureCount())
     
    # driver = ogr.GetDriverByName('Memory')
    # dataset_re = driver.CreateDataSource('')
    # layer_re = dataset_re.CreateLayer(
    #     'result',target_proj,geom_type=ogr.wkbPolygon)
    # for field in layer_jy.schema:layer_re.CreateField(field)
    
    # for feature_qh in tqdm(layer_qh):
        # geometry_qh = feature_qh.GetGeometryRef()
        # layer_jy.Clip(geometry_qh,layer_re,callback=gdal.TermProgress_nocb)
        # layer_re = layer_jy.Intersect(geometry_qh, None, 0)
        # print(layer_re.GetFeatureCount())
        # for feature in tqdm(layer_re): 
            # print(feature.items())
        # print(layer_jy.GetFeatureCount())
    
    

    

    
    # args = []
    # for feature in tqdm(layer_qh):
        # geometry = feature.GetGeometryRef()
        # wkt = geometry.ExportToWkt()
        # geometry = copy.deepcopy(geometry)
        # dataset = copy.deepcopy(dataset_jy)
        # args.append([feature_dict,geometry])
 
    # with ThreadPool(5)as pool:  # ThreadPool
    #     iterable = pool.imap(geom_intersection_layer_2_area, args)
    #     data = list(tqdm(
    #         iterable,"统计村域解译面积",len(args),disable=False))

    
    # for i in data:print(i)
    
    # for feature in tqdm(mem_lyr,desc="统计村域解译面积",disable=False):
    #     geometry = feature.GetGeometryRef()
    #     input_lyr.SetSpatialFilter(geometry)
    #     if input_lyr.GetFeatureCount() == 0 :
    #         input_lyr.ResetReading()
    #         continue
    #     areas = []
    #     for ifeature in input_lyr:
    #         igeometry = ifeature.geometry()
    #         if not(igeometry.IsSimple() and 
    #                igeometry.IsValid() 
    #                and not igeometry.IsEmpty()):
    #             igeometry = igeometry.MakeValid()
    #             igeometry = igeometry.Normalize()
    #         intersection = geometry.Intersection(igeometry)
    #         try:
    #             area = intersection.GetArea() #* 0.0015
    #             areas.append(area)
    #         except:
    #             print(feature.GetField("qhmc"))
                # print(intersection.IsSimple(),intersection.IsValid())
                # # print(ifeature.GetFID(),ifeature.GetField("OBJECTID"))
    #     input_lyr.ResetReading()
    #     record = {
    #             'mj': sum(areas) * 0.0015,
    #             'qhdm':feature.GetField('qhdm'),
    #             'qhmc':feature.GetField('qhmc'),
    #         }
    #     data.append(record)
    # input_ds.Destroy()
    # mem_ds.Destroy()
    # dbf_path = create_shp(res_shp_path,data,nf,zwdm,zwmc,jcsj,rspdid)
    # os.remove(dbf_path.replace(".dbf",".shp"))
    # os.remove(dbf_path.replace(".dbf",".shx"))
    return 


def zwcy_2_tb_(
    src_shp_path,
    dst_shp_path,
    res_shp_path,
    nf,
    zwdm,
    zwmc,
    jcsj,
    rspdid
):  
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    gdal.SetConfigOption('SHAPE_ENCODING', 'GBK') # 解决 SHAPE 文件的属性值
    input_ds = ogr.Open(src_shp_path)
    input_lyr = input_ds.GetLayer()
    target_proj = input_lyr.GetSpatialRef()
    mem_ds = transform_srs(dst_shp_path,target_proj)
    mem_lyr = mem_ds.GetLayer()
    data = []
    for feature in tqdm(mem_lyr,desc="统计村域解译面积",disable=False):
        geometry = feature.GetGeometryRef()
        input_lyr.SetSpatialFilter(geometry)
        if input_lyr.GetFeatureCount() == 0 :
            input_lyr.ResetReading()
            continue
        areas = []
        for ifeature in input_lyr:
            igeometry = ifeature.geometry()
            if not(igeometry.IsSimple() and 
                   igeometry.IsValid() 
                   and not igeometry.IsEmpty()):
                igeometry = igeometry.MakeValid()
                igeometry = igeometry.Normalize()
            intersection = geometry.Intersection(igeometry)
            try:
                area = intersection.GetArea() #* 0.0015
                areas.append(area)
            except:
                print(feature.GetField("qhmc"))
                # print(intersection.IsSimple(),intersection.IsValid())
                # # print(ifeature.GetFID(),ifeature.GetField("OBJECTID"))
        input_lyr.ResetReading()
        record = {
                'mj': sum(areas) * 0.0015,
                'qhdm':feature.GetField('qhdm'),
                'qhmc':feature.GetField('qhmc'),
            }
        data.append(record)
    input_ds.Destroy()
    mem_ds.Destroy()
    dbf_path = create_shp(res_shp_path,data,nf,zwdm,zwmc,jcsj,rspdid)
    os.remove(dbf_path.replace(".dbf",".shp"))
    os.remove(dbf_path.replace(".dbf",".shx"))
    return dbf_path


def raster_srs_transform(input_raster_path,epsg_code,output_raster_path):
    
    dst_sr = osr.SpatialReference()
    dst_sr.ImportFromEPSG(int(epsg_code))
    
    options = gdal.WarpOptions(
        xRes=10,
        yRes=10,
        dstSRS=dst_sr, 
        # srcSRS=input_sr.ExportToWkt(),
        # outputType=gdal.GDT_Float32,  
    )
    gdal.Warp(
        output_raster_path, 
        input_raster_path,
        options=options)
    # vrt_path = tif_path.replace('.tif','.vrt')
    # options = gdal.BuildVRTOptions(
    #     resolution="highest",
    #     resampleAlg=gdalconst.GRA_NearestNeighbour,
    #     separate=True)
    # gdal.BuildVRT(vrt_path,jp2_paths,options=options)
    # creationOptions = [
        # "BIGTIFF=YES",
        # "TILED=YES",
        # "BLOCKXSIZE=1024",
        # "TFW=YES",
        # "NUM_THREADS=ALL_CPUS"
    # ]
   
    # options = gdal.TranslateOptions(
    #     format="Mem",
        # xRes=10,
        # yRes=10,
        # outputSRS=dst_sr,
        # resampleAlg=gdalconst.GRA_NearestNeighbour,
        # creationOptions=creationOptions)
    
    # mem_path = "/vsimem/output.tif"
    # dataset = gdal.Translate(
    #     mem_path,
    #     input_raster_path,
    #     options=options)
    
    # creationOptions = [
    #     "BIGTIFF=YES",
    #     "TILED=YES",
    #     "BLOCKXSIZE=1024",
    #     "TFW=YES",
    #     "NUM_THREADS=ALL_CPUS"
    # ]
    # options = gdal.TranslateOptions(
    #     format="GTiff",
        # xRes=10,
        # yRes=10,
        # resampleAlg=gdalconst.GRA_NearestNeighbour,
    #     creationOptions=creationOptions)
    # gdal.Translate(
    #     output_raster_path,dataset,options=options)
    # gdal.Warp(
    #     dstSRS=output_sr.ExportToWkt(),
    #     srcSRS=input_sr.ExportToWkt(),
    # )
    
    
    # src_ds = gdal.Open(input_raster_path,gdal.GA_ReadOnly)
    # xsize = src_ds.RasterXSize
    # ysize = src_ds.RasterYSize
    # count = src_ds.RasterCount
    # srs_trans = src_ds.GetGeoTransform()
    # projectionref = src_ds.GetProjectionRef()
    # spatialref = src_ds.GetSpatialRef()
    # band = src_ds.GetRasterBand(1)
    # nodate = band.GetNoDataValue()
    # datatype = band.DataType
    
    # print(srs_trans)
    
    # dst_sr = osr.SpatialReference()
    # dst_sr.ImportFromEPSG(int(epsg_code))
    # src_sr = osr.SpatialReference(wkt=projectionref)
    # tx = osr.CoordinateTransformation(src_sr, dst_sr)
    
    # (ulx,uly,_) = tx.TransformPoint(srs_trans[0], srs_trans[3])
    # (urx,ury,_) = tx.TransformPoint(srs_trans[0] + srs_trans[1] * xsize, srs_trans[3])
    # (llx,lly,_) = tx.TransformPoint(srs_trans[0], srs_trans[3] + srs_trans[5] * ysize)
    # (lrx,lry,_) = tx.TransformPoint(
    #     srs_trans[0] + srs_trans[1] * xsize + srs_trans[2] * ysize,
    #     srs_trans[3] + srs_trans[4] * xsize + srs_trans[5] * ysize)
    
    # min_x =min(ulx,urx,llx,lrx)
    # max_x =max(ulx,urx,llx,lrx)
    # min_y =min(uly,ury,lly,lry)
    # max_y =max(uly,ury,lly,lry)
    # dst_trans = (min_x, 10, srs_trans[2], max_y,srs_trans[4], 10)
    
    # print(dst_trans)
    
    # src_wkt = src_sr.ExportToWkt()
    # dst_wkt = dst_sr.ExportToWkt()
    
    # driver = gdal.GetDriverByName('GTiff')
    # dst_ds = driver.Create(
    #     output_raster_path, 
    #     xsize, 
    #     ysize, 
    #     count, 
    #     datatype)
    # dst_ds.SetProjection(dst_wkt)
    # dst_ds.SetGeoTransform(dst_trans)
    
    # # 重投影和重采样
    # gdal.ReprojectImage(
    #     src_ds, 
    #     dst_ds, 
    #     src_wkt, 
    #     dst_wkt, 
    #     gdal.GRA_Bilinear)
 
    # # 关闭数据集
    # dst_ds = None
    # src_ds = None
    return


def raster_2_vector(tif_path,shp_path):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    dataset = gdal.Open(tif_path,gdal.GA_ReadOnly)
    projection = dataset.GetProjection()
    srs = osr.SpatialReference(wkt=projection)
    if srs.IsGeographic():print("使用的是地理坐标系")
    epsg_code = srs.GetAttrValue('AUTHORITY', 1)
    if epsg_code != '4326':
        mem_path = '/vsimem/output.tif'
        print("转WGS84坐标系")
        dateset = gdal.Warp(
            mem_path, 
            dataset, 
            dstSRS='EPSG:4326',
            resampleAlg=gdal.GRA_Bilinear,
            callback=gdal.TermProgress_nocb)
        projection = dataset.GetProjection()
        srs = osr.SpatialReference(wkt=projection)
    band = dataset.GetRasterBand(1)
    mask_band = band.GetMaskBand()
    driver = ogr.GetDriverByName('ESRI Shapefile') # 'ESRI Shapefile'
    vector_ds = driver.CreateDataSource(shp_path)
    layer = vector_ds.CreateLayer('layer',srs,geom_type=ogr.wkbPolygon)
    field_definition = ogr.FieldDefn("class", ogr.OFTString)
    field_definition.SetWidth(6)
    layer.CreateField(field_definition)
    print("栅格转矢量")
    options = [
        '8BIT_UNSIGNED=YES', 
        'ATTRIBUTE=value'
    ]
    gdal.Polygonize(
        band, 
        mask_band, 
        layer, 
        0, 
        options=options, 
        callback=gdal.TermProgress_nocb)   
    return


def zscy_2_tb(shp_path1,shp_path2,shp_path3):
    gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'YES') # 解决中文路径
    # gdal.SetConfigOption('SHAPE_ENCODING', 'UTF8') # 解决 SHAPE 文件的属性值
    # gdal.SetConfigOption('SHAPE_ENCODING', 'GBK') # 解决 SHAPE 文件的属性值
    jy_ds = ogr.Open(shp_path1)
    jy_lyr = jy_ds.GetLayer()
    
    qh_ds = ogr.Open(shp_path2)
    qh_lyr = qh_ds.GetLayer()
    
    driver = ogr.GetDriverByName('ESRI Shapefile')
    zs_ds = driver.CreateDataSource(shp_path3)
    zs_lyr = zs_ds.CreateLayer('result')
    
    field_definition = ogr.FieldDefn("mj1", ogr.OFTReal)
    field_definition.SetPrecision(24)
    zs_lyr.CreateField(field_definition)
    
    field_definition = ogr.FieldDefn("mj2", ogr.OFTReal)
    field_definition.SetPrecision(24)
    zs_lyr.CreateField(field_definition)
    
    field_definition = ogr.FieldDefn("mj3", ogr.OFTReal)
    field_definition.SetPrecision(24)
    zs_lyr.CreateField(field_definition)
    
    layerdefn = zs_lyr.GetLayerDefn()

    fields = ["mj1","mj2","mj3"]
    values = ['1','2','3']
    
    driver = ogr.GetDriverByName('Memory') # 'ESRI Shapefile'
    mem_ds = driver.CreateDataSource('')
    
    
        
    for feature in tqdm(qh_lyr,desc="统计村域级长势面积",disable=False):
        geometry = feature.GetGeometryRef()
        jy_lyr.SetSpatialFilter(geometry)
        if jy_lyr.GetFeatureCount() == 0 : 
            jy_lyr.ResetReading()
            continue
        mem_lyr = mem_ds.CopyLayer(jy_lyr, 'tmp_lyr')
        record = {}
        for field,value in zip(fields,values):
            mem_lyr.SetAttributeFilter(f"class='{value}'")
            if mem_lyr.GetFeatureCount() == 0 : 
                record.update({field:0})
            else:
                areas = []
                for ifeature in mem_lyr:
                    igeometry = ifeature.geometry()
                    if not(igeometry.IsSimple() and 
                        igeometry.IsValid() 
                        and not igeometry.IsEmpty()):
                        igeometry = igeometry.MakeValid()
                        igeometry = igeometry.Normalize()
                    intersection = geometry.Intersection(igeometry)
                    area = intersection.GetArea() 
                    areas.append(area)
                record.update({field:sum(areas) * 0.0015})
        
    
        rfeature = ogr.Feature(layerdefn)
        for k,v in record.items():rfeature.SetField(k,v)
        zs_lyr.CreateFeature(rfeature)
        
        mem_ds.DeleteLayer(0) 
        jy_lyr.ResetReading()
 
    jy_ds.Destroy()
    qh_ds.Destroy()
    
    # dbf_path = create_shp(res_shp_path,data,nf,zwdm,zwmc,jcsj,rspdid)
    # os.remove(dbf_path.replace(".dbf",".shp"))
    # os.remove(dbf_path.replace(".dbf",".shx"))
    return


if __name__ == "__main__":
    shp_path = r"/data/fengyy/dateset/temp/0429/江西村级.shp"
    check_geometry(shp_path)
    
    vector_source_path=r"/data/fengyy/dateset/temp/2024江西心连心早稻_造.shp"
    # vector_source_path=r"/data/fengyy/dateset/temp/0429/江西早稻.shp"
    vector_method_path=r"/data/fengyy/dateset/temp/36江西_村级界线.shp"
    # vector_method_path=r"/data/fengyy/dateset/temp/0429/江西村级_32650.shp"
    vector_dir=r"/data/fengyy/dateset/temp/0429"
    layer_name="2024江西村级早稻面积统计"
    kwargs = {
        'nf':"2024",
        'zwdm':"07",
        'zwmc':"早稻",
        'jcsj':"2024-03-01",
        'rspdid':"dfasdfafasdf",
    }
    # sz = Statistics_Zzjg(
    #     vector_source_path,
    #     vector_method_path,
    #     vector_dir,
    #     layer_name,
    #     kwargs
    # )
    # sz.running()

    # dataset = ogr.Open(vector_method_path,gdal.GA_ReadOnly)
    # layer = dataset.GetLayer()
    # layer.SetAttributeFilter(f"5<=FID and FID<=10")
    # print(layer.GetFeatureCount())
    
    # vector_source=r"/data/fengyy/dateset/temp/0429/江西村级.shp"
    # vector_result=r"/data/fengyy/dateset/temp/0429"
    # layer_name="江西村级_32650"
    # epsg_code=32650
    # transform_srs(
    #     vector_source,
    #     vector_result,
    #     layer_name,
    #     projection=None,
    #     epsg_code=epsg_code,
    # )
    
    # vector_source=r"/data/fengyy/dateset/temp/36江西_村级界线.shp"
    # vector_result=r"/data/fengyy/dateset/temp/0429"
    # layer_name="江西村级"
    # vector_source=r"/data/fengyy/dateset/temp/2024江西心连心早稻_造.shp"
    # vector_result=r"/data/fengyy/dateset/temp/0429"
    # layer_name="江西早稻"
    # repair_geometry(
    #     vector_source,
    #     vector_result,
    #     layer_name,
    # )

    
    # vector_source = r"/data/fengyy/dateset/temp/江西水稻.shp"
    # vector_method = r"/data/fengyy/dateset/temp/江西村域.shp"
    # vector_result = r"/data/fengyy/dateset/temp/county" # county
    # layer_name = "3602_水稻"
    # filter_params = "citycode='3602'" # "citycode='3601'" countycode='3601'
    # layer_clip_layer(
    #     vector_source,
    #     vector_method,
    #     vector_result,
    #     layer_name,
    #     filter_params
    # )
    
    # dsttif=r"/data/fengyy/dateset/temp/L1C_50SKA_20240329-053813.SAFE_1.tif"
    # srctif=r"/data/fengyy/dateset/temp/L1C_50SKA_20240329-053813.SAFE_0.tif"
    # shppath=r"/data/fengyy/dateset/temp/cy_41.shp"
    # shp_clip_tif(dsttif,srctif,shppath)
    
    # src_shp_path=r"/data/fengyy/dateset/temp/2024江西心连心早稻_造.shp"
    # dst_shp_path=r"/data/fengyy/dateset/temp/36江西_村级界线.shp"
    # res_shp_path=r"/data/fengyy/dateset/temp/江西早稻村级.shp"
    # nf="2024"
    # zwdm="07"
    # zwmc="早稻"
    # jcsj="2024-03-01"
    # rspdid="dfasdfafasdf"
    # zwcy_2_tb(src_shp_path,dst_shp_path,res_shp_path,nf,zwdm,zwmc,jcsj,rspdid)
    
    # vector_source_path = r"/data/fengyy/dateset/temp/江西水稻.shp"
    # vector_method_path = r"/data/fengyy/dateset/temp/江西村域.shp"
    # data = shp_clip_shp(vector_source,vector_method)
    # load(vector_source_path,vector_method_path)
    
    # savedir=r"/data/fengyy/dateset/temp"
    # layername = r"江西村域水稻统计"
    # kwargs = {
    #     "nf": "2024",
    #     "zwdm": "07",
    #     "zwmc": "早稻",
    #     "jcsj": "2024-03-01",
    #     "rspdid": "dfasdfafasdf",
    # }
    # data_2_dbf(data,savedir,layername,kwargs)
    
    # shp_path = r"/data/fengyy/dateset/temp/36江西_村级界线.shp"
    # out_path = r"/data/fengyy/dateset/temp/36江西_村级界线_修复.shp"
    # dataset = repair_geometry(shp_path)
    # layer = dataset.GetLayer()
    # print(layer.GetFeatureCount())
    
    

# src_shp_path = r"/data/fengyy/dateset/temp/36江西_村级界线.shp"
# epsg_code=32650
# save_dir = r"/data/fengyy/dateset/temp"
# layer_name = "36江西_村级界线_32650"
# transform_srs_zipfile(src_shp_path,epsg_code,save_dir,layer_name)

# shp_path = r"/data/fengyy/dateset/temp/36江西_村级界线_32650/36江西_村级界线_32650.shp"
# check_geometry_vaild(shp_path)

# tif_path = r"/data/fengyy/dateset/temp/2023周口市项城市4乡镇玉米0716长势/2023周口市项城市4乡镇玉米0716长势3级.tif"
# shp_path = r"/data/fengyy/dateset/temp/2023周口市项城市4乡镇玉米0716长势/result.shp"
# raster_2_vector(tif_path,shp_path)


# shp_path1 = r"/data/fengyy/dateset/temp/2023周口市项城市4乡镇玉米0716长势/result.shp"
# shp_path2 = r"/data/fengyy/dateset/temp/cy_41.shp"
# shp_path3 = r"/data/fengyy/dateset/temp/2023周口市项城市4乡镇玉米0716长势/test.shp"
# zscy_2_tb(shp_path1,shp_path2,shp_path3)


# input_raster_path = r"/data/fengyy/dateset/temp/2023周口市项城市4乡镇玉米0716长势/2023周口市项城市4乡镇玉米0716长势3级.tif"
# epsg_code = 32650
# output_raster_path = r"/data/fengyy/dateset/temp/2023周口市项城市4乡镇玉米0716长势/test.tif"
# output_raster_path = r"/vsimem/output.tif"
# raster_srs_transform(input_raster_path,epsg_code,output_raster_path)




# src_shp_path=r"/data/fengyy/dateset/temp/2024三亚市早稻20240115.shp"
# epsg_code = "32649"
# save_dir=r"/data/fengyy/dateset/temp"
# layer_name = "2024三亚市早稻20240115_49rn_result"
# transform_srs_zipfile(src_shp_path,epsg_code,save_dir,layer_name)
